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ABSTRACT 

When  a  marine  vessel  is  used  as  a  piafform  for  tracking  space  vehicles,  the  accuracy 
af  the  tracking  is  critically  limited  by  errors  in  the  ship's  estimated  position  and  ori¬ 
entation.  These  reference  data  are  usually  estimated  by  combining  the  output  of  an 
inertial  navigation  system  with  additional  external  position,  velocity  on^/or  orien¬ 
tation  measurements.  To  meet  stringent  accuracy  requirements  for  this  estimation, 
optimum  statistical  computation  procedures  for  combining  the  data  are  required. 

This  report  develops  a  discrete-time,  linear  model  for  the  error  propagation  in  an  in¬ 
ertial  navigation  system.  This  model  mokes  it  possible  to  derive  on  optimum  recursive 
data  processing  procedure  for  combining  the  inertial  system  output  with  the  external 
fixes.  A  recursive  formula  for  the  minimum-meon-squore  error  in  estimated  position, 
velocity  and  orientation  is  also  obtained.  The  various  formulas  ore  directly  appli¬ 
cable  to  broad  classes  of  inertiol  navigation  systems  and  external  measurements. 
A  general  discussion  is  provided  on  the  computational  problems  associated  with  both 
actual  data  processing  and  the  performance  of  system-error  analyses. 
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I.  INTRODUCTION 

The  problem  that  motivated  the  analysis  contained  in  this  report  is  connected  with  the  use 
of  marine  vessels  as  platforms  for  devices  used  to  track  missiles,  satellites  and  interplanetary 
probes.  One  of  the  major  difficulties  associated  with  these  marine -stationed  tracking  platforms 
is  the  degradation  of  the  usefulness  of  the  tracking  data  caused  by  uncertainties  as  to  the  ship's 
latitude,  longitude  and  orientation  with  respect  to  north  and  vertical.  Techniques  that  can  be 
used  to  determine  these  reference  data  for  a  marine  vessel  can  be  subdivided  into  two  basic 
classes.  In  the  first  class  is  the  use  of  direct  measurements  such  as  star  fixes,  tracking  of 
satellites,  shore-based  electromagnetic  navigation  aids  and  underwater  sonar  beacons.  In  the 
second  class,  consisting  of  direct  navigational  procedures,  the  technique  of  prime  interest  is  the 
inertial  navigation  system.^  These  two  methods,  external  measurements  and  inertial  navigation 
systems,  have  different  performance  characteristics.  An  inertial  navigation  system  is  capable 
of  high  accuracy  for  short  periods  of  time,  but  the  errors  can  build  up  to  an  unacceptable  level. 
External  measurements  provide  information  only  at  the  time  of  measurement,  and  any  single 
measurement  might  not  be  accurate  enough  to  meet  requirements  or  might  not  give  enough  infor¬ 
mation  to  completely  specify  the  desired  aspects  of  the  vessel's  position,  velocity  and  orienta¬ 
tion.  Because  of  economic  as  well  as  physical  limitations,  the  required  accuracies  are  not  al¬ 
ways  achievable  by  merely  improving  the  inertial  navigator  or  the  basic  measurement  systems. 

When  combined  in  an  efficient  manner,  however,  the  two  techniques  complement  each  other 
and  can  possibly  provide  a  resulting  capability  far  superior  to  either  single  procedure.  The  in¬ 
ertial  navigator  allows  the  use  of  incomplete  external  measurements  taken  at  different  times, 
and  statistical  smoothing  of  the  random  errors  in  the  external  measurements  can  be  achieved. 

In  turn,  the  external  measurements  help  control  the  error  build-up  inherent  in  the  inertial  nav¬ 
igation  system.  Since  extensive  computation  facilities  for  data  processing  are  usually  available 
to  the  marine -stationed  tracking  platform,  one  can  use  optimum  data  processing  techniques  for 
combining  external  measurements  with  the  data  obtained  from  an  inertial  navigation  system  to 
provide  a  best  estimate  of  a  ship's  position  and  orientation.  Knowledge  of  such  optimum  compu¬ 
tation  procedures  also  provides  answers  to  such  related  operational  questions  as:  What  external 
measurements  are  to  be  taken,  at  what  time,  and  in  what  combinations  to  provide  the  most  effec¬ 
tive  results?  When  can  inexpensive  external  measuring  devices  be  combined  with  an  inexpensive 
inertial  navigation  system  to  provide  performance  characteristics  competitive  with  far  more 

tTha  inertial  navigation  system  could  be  considered  os  using  external  measurements  in  that  it  relies  on  gravita¬ 
tional  forces.  However,  the  dichotomy  is  useful  for  the  sake  of  exposition. 
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expensive  inertial  navigation  systems?  If  underwater  sonar  beacons  are  being  employed,  how 
should  they  be  placed  and  in  what  pattern  about  them  should  the  marine  vessel  sail? 

This  report  develops  a  theory  and  derives  explicit  formulas  that  provide  answers  to  the  above 
questions.  The  results  are  general  in  nature  and  thus  are  applicable  to  problems  beyond  those 
that  were  of  interest  initially.  A  discrete-time  model  for  the  error  propagation  of  a  marine  in¬ 
ertial  navigation  system  is  developed  and,  on  the  basis  of  this  model,  an  optimum  method  for 
combining  the  external  measurements  with  the  inertial  navigation  system  is  derived.  The  opti¬ 
mum  technique  provides  a  reasonable  and  workable  method  of  data  processing,  but  even  if  re¬ 
strictions  prevent  its  explicit  use,  the  theory  provides  a  valuable  standard  of  comparison,  as 
well  as  a  good  starting  point  for  modifications.  The  derivation  of  the  data  processing  procedure 
also  includes  derivation  of  formulas  that  determine  how  the  various  errors  in  the  inertial  naviga¬ 
tion  system  and  external  measurements  are  propagated  through  the  optimum  data  processing  tech¬ 
nique.  These  formulas  thus  provide  a  very  powerful  tool  for  comparing  and  evaluating  the  effec¬ 
tiveness  of  different  external  measurement  and  inertial  navigation  system  configurations. 

Explicit  formulas  are  provided,  but,  because  of  the  inherent  complexity  of  the  error  behav¬ 
ior  of  inertial  navigation  systems,  closed-form  solutions  are  usually  not  possible.  The  formulas, 
however,  are  directly  solvable  on  a  digital  computer  (some  analog  equipment  could  also  be  used 
if  desired)  by  the  use  of  straightforward  techniques.  For  example,  the  formula  governing  the  er¬ 
ror  propagation  in  time  is  a  nonlinear  recursive  equation  whose  solution  can  be  calculated  di¬ 
rectly  in  time  from  known  initial  conditions.  Since  the  necessary  computer  programs  have  not 
been  written,  no  explicit  numerical  results  are  included  in  this  report. 

The  material  in  this  report  assumes  familiarity  with  inertial  navigation  systems;  therefore, 
the  treatment  of  the  systems  themselves,  although  mathematically  complete,  is  very  abbreviated. 
The  development  of  the  theory  for  the  optimum  techniques  is  more  expository,  but  assumes  a 
basic  knowledge  of  random  variables  and  stochastic  processes. 

n.  GENERAL  DISCUSSION 

Because  of  the  technical  nature  of  the  theory  and  formulas  to  be  developed,  it  is  appropriate 
to  preface  the  actual  derivations  with  a  general  discussion  of  the  basic  assumptions  and  methods 
of  analysis  to  be  used.  The  remainder  of  the  report  (Secs.  Ill,  IV  and  V)  will  be  of  interest  pri¬ 
marily  to  those  readers  who  wish  to  employ  the  general  theory  or  use  the  explicit  formulas  to 
solve  particular  problems. 

The  model  of  the  inertial  navigation  system  to  be  considered  includes: 

(a)  The  possible  incorporation  of  a  velocity  log  for  measuring  the  vessel's 
velocity  and  the  effect  of  errors  in  these  measurements. 

(b)  The  effect  of  errors  due  to^  gyro  drift  and  accelerometer  bias. 

(c)  Incorporation  of  an  extremely  broad  class  of  external  measurements, 

(d)  The  propagation  of  all  seven  components  of  error,  two  location  coordinates 
and  their  associated  velocities  and  three  attitude  angles.  (The  term  "posi¬ 
tion  estimation,"  as  used  throughout  the  remainder  of  this  report,  actually 
refers  to  all  seven  quantities,  not  just  the  location  of  the  ship  on  the  sur¬ 
face  of  the  earth.) 

Section  III  develops  a  discrete-time  model  of  the  error  propagation  in  the  inertial  navigator. 

In  deriving  the  error  model,  it  is  assumed  that  the  other  accelerations  acting  on  the  platform  are 
small  compared  with  the  acceleration  of  the  earth's  gravitational  field  and  that  the  velocity  of  the 
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ship  with  respect  to  the  earth's  surface  is  small  compared  with  the  velocity  of  the  earth's  sur¬ 
face  with  respect  to  some  inertial  reference.  We  assume  that  the  velocity-log  errors  and 
accelerometer-bias  errors  (additive  errors  in  the  accelerometer  output)  can  be  represented  by 
continuous,  stationary,  zero  mean,  Markov  random  processes.  Further,  it  is  assumed  that  the 
values  of  the  associated  correlation  functions  for  time  shifts  larger  than  the  time  between  fix  ob¬ 
servations  will  be  very  small  compared  with  the  values  for  zero  time  shift.  Stationary  first-order 
Markov  processes  are  assumed  for  the  gyro  drifts,  but  here  it  is  not  assumed  that  the  correla¬ 
tion  times  are  small  with  respect  to  the  time  between  fix  observations.  The  linearized  model 
for  the  error  propagation  is  of  the  form 

x(t  +  T)  =  *(t  +  T)x(t)  +  w(t)  ,  (1) 

in  which  T  is  the  time  between  external  measurements.  The  column  vector  x(t)  is  a 
ten-dimensional  vector,  the  first  seven  entries  of  which  are  the  seven  error  variables  while 
the  last  three  are  variables  that  are  the  discrete-time  equivalent  of  the  gyro-drift  inputs  (equiv¬ 
alent  in  the  sense  that  they  produce  the  same  response  at  all  sample  times  as  the  original  con¬ 
tinuous  gyro-drift  inputs) .  The  column  vector  w(t)  is  also  a  ten-dimensional  vector,  the  first 
seven  entries  of  which  are  inputs  that  directly  represent  the  effects  of  accelerometer  bias  and 
velocity-log  errors.  The  last  three  entries  of  w(t)  are  inputs  that,  when  passed  through  a  first- 
order  system,  result  in  the  equivalent  gyro-drift  inputs  (the  last  three  entries  of  x).  The  10  X  10 

matrix  ♦  thus  includes  both  the  transition  matrix  of  the  error-propagation  system  and  the  tran- 
st 

sition  matrix  of  the  1  -order  system  used  to  generate  the  three  equivalent  gyro-drift  inputs. 

The  input  w(t)  is  white  in  the  sense  thatt 

®^^t-t+kT^  ~  ^  ^  ■  integer;  k  ^  0  , 

where  the  superscript  T  denotes  transposition.  The  matrix  ♦  is  a  constant  for  picket  ship  op¬ 
eration  (operation  confined  to  a  distance  of  several  miles  from  a  fixed  point)  and  varies  with  lo¬ 
cation  for  normal  marine  operation. 

Section  IV  treats  the  estimation  problem.  For  the  external  measurements,  we  assume  a 
rather  general  form 

r(t)  =  F[q(t)]  -  v(t)  , 

in  which  the  vector  x  is  the  observation  taken,  q(t)  is  the  seven-vector  giving  the  true  quantities 
associated  with  the  inertial  navigator  operation,  and  v(t)  is  the  measurement  error.  This  error 
is  assumed  to  have  zero  mean  and  to  be  white  in  the  sense  that 

E{v^v^^j^rp}  =0  k  =  integer;  k  ^0 

The  output  of  the  inertial  navigator  is  jr(t)  =  q(t)  +  x(t).  The  variables  used  by  our  over-all  sys¬ 
tem  to  estimate  x^  are 

z(t)  =  F(^(t)]  -  r(T) 

=  H(t)x(t)  +  y(T)  T  =  t  —  kT,  k  =  0,  1,  Z,  .  .  .  ,  (2) 

tThe  symbol  x(t)  denotes  a  sample  function  of  a  random  process  (a  function  of  time  that  is  the  result  of  a  single 
realization,  or  experiment)  while  the  symbol  x^  denotes  a  random  variable  (on  observation  made  at  a  single  fixed 
time  that  is  a  function  of  the  realization,  or  experiment  observed).  Thus,  a  time  average  is  defined  only  for  x(t) 
and  a  statistical  average  only  for  x^. 
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in  which  we  have  expanded  F[^(t)]  in  a  Taylor's  series  about  ^{r)  =  q(T);  that  is, 
9F,[s(t)]| 
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y(T)  =  q(T) 

Thus,  the  external  measurements  are  also  linearized. 

The  use  of  a  linearized  model  such  as  that  represented  by  Eqs.  (1)  and  (2)  is  justified  in  most 
problems  of  interest,  since  we  are  primarily  concerned  with  highly  accurate  systems  having  small 
errors.  Of  course,  in  any  particular  application,  one  must  always  determine  whether  or  not  the 
linearized  model  is  a  valid  representation  of  the  important  aspects  of  the  physical  situation. 

The  estimation  scheme  developed  in  Sec.  IV  is  linear;  that  is,  the  optimum  estimate  is  ex¬ 
pressed  as 


£(t)  =  2  Ck^(t  -  k)  . 

in  which  the  sum  extends  over  all  measurements.  The  reasons  for  using  a  linear  estimation 
scheme  are  twofold: 

(a)  The  solution  to  the  linear  estimation  problem  is  tractable. 

(b)  It  is  unrealistic  to  assume  that  we  would  have  good  estimates  of  more  than 
second-order  statistics;  using  only  second-order  statistics,  the  linear  es¬ 
timate  is  optimum. t  Moreover,  it  seems  reasonable  to  assume  that  the 
sources  generating  navigator  errors  would  be  Gaussian.  Since  the  error- 
propagation  system  is  linear  to  a  good  approximation,  the  processes  that 
we  are  trying  to  estimate  may  be  assumed  Gaussian  to  a  good  approxima¬ 
tion.  Under  these  conditions,  the  optimum  mean-square  estimator  is  a 
linear  estimator. 
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The  solution  to  the  estimation  problem  is  similar  to  Kalman's  recursive  formulation.  The  re¬ 
sults  are; 

(a)  A  recursive  formula  that  provides  the  optimum  data  processing  procedure. 

(b)  A  recursive  relation  for  determining  the  covariance  matrix  of  the  error 
in  estimating  the  position  of  the  marine  vessel. 

The  recursive  formulations  given  in  Sec. IV  are  nonlinear  and,  for  almost  all  situations  of 
interest,  require  machine  solution.  Sec.V  discusses  the  computational  steps  necessary  to  carry 
out  a  machine  solution. 

One  final  point:  the  formulas  to  follow  are  complex  because  of  the  generality  of  the  assumed 
model.  In  any  explicit  application,  extensive  simplifications  may  be  possible.  Similarly,  the 
analysis  itself  can  be  extended  in  various  directions  by  using  the  same  basic  techniques  employed 
in  this  report.  To  cite  just  one  example,  the  theory  is  directly  extendible  to  external-fix  meas¬ 
urement  errors  that  are  correlated  in  time  and  that  have  a  nonzero  bias  component. 

in.  MODEL  OF  THE  NAVIGATOR  ERROR  PROPAGATION 

In  this  section  we  proceed  directly  to  the  problem  of  deriving  a  model  for  the  error  propaga¬ 
tion  in  a  marine  inertial  navigation  system.  The  reader  unfamiliar  with  inertial  navigation  sys¬ 
tems  may  wish  to  consult  the  literature,  since  certain  of  our  analyses  are  succinct  rather  than 


t  We  consider  the  optimum  or  best  estimate  of  position  to  be  the  one  that  simultaneously  minimizes  the  voriance 
of  each  component  of  position.  It  is  shown  in  Sec.  IV  that  this  estimate  also  yields  an  ellipsoid  of  concentration^ 
that  falls  inside  the  ellipsoid  of  concentration  of  any  other  estimate. 
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expository.  There  is  at  least  one  textbook  which  covers  the  general  subject,  although  the  author 
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found  an  unpublished  set  of  notes  to  be  somewhat  more  informative.  The  idea  behind  an  inertial 
navigation  system  is  quite  basic.  A  gyro-stabilized  platform  is  used  to  orient  two  accelerometers. 
These  sense  components  of  acceleration  in  the  east  and  north  directions,  and  these  accelerations 
are  integrated  twice  by  a  computer  to  obtain  position.  The  computer  also  computes  the  torques 
necessary  to  keep  the  platform  orientation  north  and  east  as  the  earth  rotates  and  the  ship  moves. 

Many  marine  inertial  navigation  systems  also  include  velocity- log  feedback  in  which  estimates 
of  the  vessel's  velocity,  obtained  from  an  electromagnetic  velocity  meter,  are  fed  into  the  system 
to  provide  damping  of  certain  errors.  Although  such  velocity  information  could  be  considered  ex-- 
ternal  measurements,  the  analyses  to  follow  include  velocity-log  feedback  as  an  intrinsic  part  of 
the  inertial  navigation  system.  It  will  be  seen,  however,  that  this  approach  greatly  complicates 
the  equations.  The  possibilities  of  treating  the  velocity- log  data  as  external  measurements  are 
briefly  discussed  in  Sec.  V. 


A.  Definitions  of  Symbols  and  Coordinate  Systems 

Throughout  this  report,  vector  quantities  are  indicated  by  underlining.  The  conventional 
latitude-longitude  system  is  used;  that  is, 

X  =  degrees  north  (+)  latitude 

and 


A  =  degrees  east  (+)  longitude 

are  the  variables  describing  the  ship's  position.  The  operation  of  the  navigation  system  will  be 
based  on  a  set  of  coordinates  whose  origin  is  at  the  center  of  the  earth  and  whose  z-axis  passes 
through  the  center  of  the  ship.  These  coordinates  will  be  oriented  with  the  x-axis  pointed  north 
and  the  y-axis  pointed  west.  In  the  remainder  of  this  report,  this  coordinate  system  will  be  re¬ 
ferred  to  simply  as  the  system  coordinates.  We  will  also  need  to  distinguish  between  three  such 
sets  of  system  coordinates; 

(1)  X,  y,  z:  The  set  of  system  coordinates  based  on  the  true  position  of  the  ship. 

(2)  Xq,  y^,  z^\  The  set  of  system  coordinates  whose  orientation  is  based  on  the 
position  of  the  ship  calculated  by  the  navigation  system  computer. 

(3)  xp,  yp,  zp:  The  set  of  system  coordinates  aligned  with  the  orientation  of  the 
stable  platform,  i.e.,  x^  parallel  to  x  accelerometer. 

The  difference  between  the  last  two  sets  of  coordinates  is  due  to  the  curvature  of  the  earth  and 
the  existence  of  error  in  the  navigator-calculated  position. 

Note  that  these  three  coordinate  systems  are  all  related  to  one  another  by  rotational  transfor¬ 
mations.  A  rotational  transformation  is  specified  by  three  angles.  If  all  these  angles  are  small 
and  sines  are  approximated  by  the  angle  and  cosines  by  one,  then  the  rotational  transformation 
can  be  (approximately)  specified  by  a  linear  transformation  in  which  these  angles  appear  as  coeffi¬ 
cients.  We  assume  that  our  error  angles  are  small.  Thus,  the  disparity  among  the  three  coor¬ 
dinate  systems  can  be  described  by  the  vectors; 

£  =  vector  error  angle  between  the  true  axes  and  platform  axes, 

86  =  vector  error  angle  between  the  true  axes  and  computer  axes, 

^  =  vector  error  angle  between  computer  axes  and  platform  axes, 

<p  =  ip  +  86. 
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Each  of  these  vector  quantities  has  three  components;  e.g.,  j/i^has  the  components  ili^,  iji^, 
representing  rotations  about  the  x,  y,  and  z  axes^  respectively. 

We  also  introduce  the  following  quantities: 

A  =  acceleration  sensed  by  the  accelerometer. 

R  =  position  vector  of  the  vehicle  as  measured  in  a  coordinate  system  aligned 
with  X,  y  and  z  whose  origin  is  at  the  center  of  the  earth,  i.e.,  the  system 
coordinates.  Note  that  x  =  y  =  0,  z=R,  R  =  earth's  radius  =  |r|  . 

R  =  velocity  vector  of  vehicle  in  the  system  coordinates.  Note  that  R  =  0. 
gj^(R)  =  mass  attraction  of  the  earth  at  position  R. 

w  =  angular  velocity  of  system  coordinates  relative  to  an  inertial  set  of  coor¬ 
dinates,  that  is,  a  set  of  coordinates  stationary  with  respect  to  the  fixed 
stars. 

n  =  angular  velocity  of  earth. 

n  =  frequency  of  earth's  angular  rotation  about  the  North  Pole  =  |  111 . 

=  error  in  latitude  =  X-  —  Xa/.»„ai.  in  which  X  is  the  computed  latitude.  Note 
that  Ax  =  RAX,  Ay=-R  cSsXAA. 

AX  =  derivative  of  AX  with  respect  to  time. 

=  error  in  x-component  of  ship's  velocity  log  [an  electromagnetic  (actually 
magnetohydrodynamic)  device  for  measuring  the  ship's  velocity]. 

K  =  feedback  from  velocity  log  (units  of  inverse  seconds) . 

6A^  =  bias  error  in  the  x-component  of  the  accelerometer  output;  that  is,  is 
a  random  process  which  constitutes  the  additive  error  in  the  x-accelerom- 
eter  output. 

=  x-component  of  gyro  drift. 

Similar  definitions  apply  to  BV  ,  6A  ,  e  and  e  . 

y  y  y  z 

In  the  next  two  sections  we  derive  the  differential  equations  governing  the  error  propagation 

in  the  navigator.  These  equations  can  be  divided  into  two  groups.  The  first  group  is  the  set  of 

differential  equations  for  AX  and  AA  in  which  *  ,  ip  and  ^  appear  as  driving  terms.  The  second 

X  y  z 

group  of  equations  describes  the  propagation  of  the  terms  ip  ,  ip  and  ip  . 

X  y  z 

B.  Equations  Governing  the  Propagation  of  AX  and  A  A 

Recalling  the  method  of  operation  of  the  inertial  navigator,  we  see  that  the  fundamental  equa¬ 
tion  describing  the  system  is 

• 

in  which  I  indicates  that  the  differentiation  is  carried  out  with  respect  to  an  inertial  set  of  coor¬ 
dinates.  But,  using  the  general  expression  relating  differentiation  in  two  coordinate  systems 
whose  relative  motion  is  only  rotation,  we  obtain 

(^)i  =  R  + «  .  W 

in  which  w  x  R  is  the  cross  product  of  the  vectors  u  and  R.  Also, 
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Substituting  Eq.  (4)  into  Eq.  (5),  we  obtain 


^  Sya  Coord  ^  (R +i-  X  R) 

=  R  +  20)  X  R  +  u  X  (u  X  R)  +  w  X  R  .  (6) 

Substituting  Eq.(6)  into  Eq.  (3)  and  recalling  that  R  =  0,  we  obtain 

A  = -gj^(R)  +  0)  X  (w  X  R)  +  o)  X  R  .  (7) 

Now, 

x  =  y=0,  z  =  R 

Further, 

0)^  =  (41+  A)  cos  X,  o)y  =  X,  o»^  =  (II  +  A)  sinX 
Thus  Eq.  (7)  becomes 

=  RX  +  R(I1  +  A)^  COS.X  sinX  —  g^(R) 

Ay  =  —  R  cos  XA  +  2R  (12  +  A)  X  sin  X— ^(R)  .  (8) 

The  computation  scheme  used  by  the  computer  is  to  calculate  the  terms 

R(I2  +  A)^  cos  X  sinX,  8jj(R)»  2R(I2  +  A)  X  sin  X,  ^^R) 

and  subtract  them  from  the  output  of  the  accelerometer,  thus  obtaining  RX  and  R  cosXA.  These 
are  integrated  twice  to  calculate  position.  The  velocity-log  data  is  incorporated  into  a  simple 
feedback  loop  to  damp  certain  errors.  The  block  diagram  of  the  navigation  system  is  shown 
in  Fig.  1. 

Writing  out  the  equations  for  the  block  diagram,  we  have  the  equations  which  govern  the  er¬ 
ror  propagation: 

R(X  +  AX)  =  RX  +  R  ((12  +  A)^  sinX  cos X- (12  +  A  +  AA)  sin  (X  +  AX)  •  cos  (X  +  AX)] 

-KRAX  +  K6V^-(^()yA^-(()^Ay)-[gJR)-g^(R  + AR)]  +  6A^  ,  (9) 

-R  cos(X  +  AX)  (A  +  AA)  =-RA  cos  X  +  2R  [(12  +  A)  X  sin  X  -  (12  +  A  +  AA)  (X  +  AX) 

X  sin  (X  +  AX)]  +  KR  cos  (X  +  AX)  (A  +  AA)  -  KR  cos  X  (A  +  (6Vy)/R] 

-  *x\'>  *  -  «y<^  +  ^)  +  6Ay]  (10) 

We  assume  that  AX  «  1  and  AA  «  1;  hence, 

sin  (X  +  AX)  »  sin  X  +  AX  cos  X, 
cos  (X  +  AX)  w  cos  X  —  AX  sin  X. 

Now  we  assume  that  the  magnitude  of  the  earth's  gravitational  attraction  at  sea  level  is  known 
for  the  region  of  operation.  If  this  is  so,  then  direction  is  the  only  unknown  quantity,  and 

g(R)  -  £(R  +  ^)  = -g^ 
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Block  diagram  of  latitude-longitude  inertial  navigation  system 


Hence, 

-gx<®  +  +  ^)  =  .  (11) 

-gy(B)  +  gy(R  +  ^)  =  -gAA  .  (12) 

It  should  be  noted  that  we  have  assumed  a  spherical  earth.  Although  the  earth's  ellipticlty 
cannot  be  ignored  in  the  actual  operation  of  the  navigation  system,  the  manner  in  which  it  affects 
error  propagation  is  secondary  and  will  be  ignored.  Similarly,  for  marine  operation,  we  can 
make  the  further  assumptions: 


KA  and  RA  are  negligible  with  respect  to  g 


We  also  note  that 


_2  1.34 


cycles/ sec  ,  cycles/ sec  ; 

"  10 ' 


hence. 


We  also  assume  that  AX  and  A  A  are  of  the  same  order  of  magnitude,  and  that  and  are 

all  of  the  same  magnitude.  If  we  apply  this  assumption  and  the  approximations  (11)  through  (14) 
to  Eqs.(9)  and  (10),  we  obtain,  after  dropping  second-order  terms, 

•  •  *  *  2  K  tf 

AX  +  2C^a)2AA  +  KAX  +  w^AX  =  ^ '^‘y  ' 

^  *  '  2  K  ft 

AA  -  2  AX  +  KAA  +  w^AA  =  -  6V  -  ,  (16) 

1  1^1 


in  which 


=  Schuler  Frequency  =  , 


"2  =  <^2”  '  "i=Ci«  * 

Cj  =  cos  X  , 


=  sin  X 

If  we  are  considering  picket  ship  operation,  then  and  are  constant  and  our  error  propaga¬ 
tion  is  determined  by  a  set  of  linear  differential  equations  with  Constant  coefficients.  For  usual 
marine  use  we  have  linear  differential  equations  with  coefficients  that  vary  with  position. 

The  quantities  ip  and  ip  act  as  drives  to  Eqs.  (15)  and  (16).  We  must  next  write  the  equations 
X  y 

that  determine  the  propagation  of  the  error  angle 


C.  Equations  CSoverning  the  Propagation  of  >p 

The  derivations  of  this  section  follow  closely  those  of  Kachickas.^  To  start,  let  us  recall 
that  we  have  defined  the  three  sets  of  coordinate  axes; 
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* 


x,y,  z: 


true  position  coordinate  axes 


X  ,  y  ,  z  :  platform  coordinate  axes 
P  P  P 

*c'^c'^c' 


All  three  coordinate  axes  have  their  origins  at  the  true  position  of  the  ship,  but  the  three  axes 
are  not  perfectly  aligned  with  one  another. 

Let 


X  denote  an  arbitrary  vector  in  the  x,  y,  z  system  , 
X  denote  the  same  vector  in  the  x  ,  y  ,  z  system  , 

-p  p’^'p'  p 

x^  denote  the  same  vector  in  the  x^,  y^,  system 

Then,  to  a  first-order  approximation, 

X  =  *x  and  x  =  0x  , 

-p  -  -c  -  ' 

in  which 


■  1 

0 

y 

•  1 

60 

z 

-60 

y 

*  = 

-‘f’z 

1 

0 

,  e  = 

-60^ 

1 

60 

X 

<p 

y 

-^x 

1 

60 

y 

-60 

X 

1 

We  have  previously  defined 


0 

[60 

X 

f.  ^ 

0 

y 

,  60  = 

60 

y 

Now,  by  definition, 

Wp  =  w  +  ^  ■  (17) 

The  platform  velocity  may  also  be  expressed  in  terms  of  the  gyro  drift  c  and  the  angular 
velocity  which,  according  to  the  computer  calculations  must  be  applied  to  the  platform  to 
maintain  its  north,  east,  up  orientation.  In  terms  of  the  true  coordinates,  this  angular  velocity 
is 

e“^a)  +  « 

—  c  — 

Hence, 
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(18) 


.1 

0)  =  #[6  Ct)  +  €] 

— p  ^  — c 

=  (*e"^)  Me  +  • 

By  using  the  given  expressions  for  *  and  0  and  noting  that,  to  a  first-order  approximation, 

1  -60  60 
z  y 


0 


-1 


60 

z 


-60  60 
y  X 


-60. 


we  obtain 


til  =0)  +  1(0  —  60)  Xw  +  €  +  0  x  e 

-p  — c  '2.  -  _c  -  -2-  - 


Now,  again  by  definition. 

Substituting  Eq.(20)  into  Eq.(i9)  and  ignoring  second-order  terms,  we  obtain 
ojp  =  w  +  ^  +  (^  —  60)  X  € 

By  equating  this  expression  to  Eq.  (17)  we  obtain 

<p  —  60  +  {0—  60)  =  €  , 

or 

>l>  +  U  X  f  =  € 


(19) 


(20) 


(21) 


(22) 


Expanding  Eq.  (22)  and  substituting  in  the  expressions  for  a>„,  and  w_,  given  following 

X  y  z  . 

Eq.  (7),  we  obtain 

(11  +  A)  C2(l'y  =  Cx  , 

(fy  +  C2(n  +  A)  C^(n  +  k)  , 

+ +  k)  . 

If  we  ignore  A  and  X  with  respect  to  il  and  assume  that  ^  ^  and  ^  are  of  the  same  magnitude, 

X  y  z 

the  above  equations  become 


w  —  W^lll  =  « 

X  2’^y  X 


V^"2^x-"l^z=  V 


»  +  u^tli  =  e 
z  2  y  X 


(23) 

(24) 

(25) 


in  which 


«2  =  <==2" 
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Fig.  2.  Flow  graph  expressing  the  dynamics  of  the  error  propagation. 


l>-tt-4lll| 

(s‘+n')*«,Ms4+KS+i.,‘) 

- 5 - 

's(s*+o*)* 

C,D 

-zc,i.i,<.>,>s(s>+a>)> 

- p - 

«,'(S'+0')'(S*  +  KS+w,') 

- 6 - 

D»  [s*  4-n*]*[(8' ♦  KS+»,*I*-4«»,'S*] 


Fig.  3.  Reduced  flow  graph. 


D.  Complete  Error  Propagation  Dynamics 

The  five  equations  (15),  (16),  (23),  (24)  and  (25)  completely  describe  the  propagation  of  errors 
in  the  navigation  system.  Our  first  objective  will  be  to  obtain  the  transfer  functions  and  impulse 
responses  relating  AX  and  AA  to  the  various  inputs.  This  can  be  accomplished  by  first  drawing 
the  flow  graph  expressing  the  five  equations: 


AX  +  2C,(i>,AA  +  KAX  +  a)  AX  =  ^ 

1  fa  S  i\ 

y 

tEq.(15)) 

w,  . 

AA  -  2  AX  +  KAA  +  0)  AA  =  - 

s 

„  6A 

[Eq.(16)] 

'^x-“2^=  "x 

[Eq.(23)] 

^y  +  «2^x-“l^z=^  ‘y 

[Eq.(24)] 

'^z  +  “2’^y="x  • 

[Eq.(25)] 

This  flow  graph  is  shown  in  its  original  form  in  Fig.  2  and  in  reduced  form  in  Fig.  3;  the  details 
of  the  reduction  are  omitted. 

In  order  to  proceed  further,  we  must  factor  the  denominator  term  D,  which  appears  in  the 
transfer  functions  in  Fig.  3; 

D  =  (S^  +  n^)^  ((S^  +  KS  +  ojg)^  -  4a»2S^] 

ai(S^  +  n^)^  [(S^  +  KS  +  1.0821  (S^  +  KS  +  0.9241  •  (26) 

s  s 

Now,  from  this  we  could  obtain  the  necessary  transfer  functions  directly.  The  resulting  expres¬ 
sions,  however,  would  be  extremely  unwieldy.  For  this  reason  we  consider  approximating  D  by 

D  «  (S^  +  n^)^  (S^  +  KS  +  Wg)^  .  (27) 

In  order  to  estimate  the  error  involved,  consider  the  time- response  terms  that  would  result  from 
each  expression.  From  Eq.  (26)  we  would  obtain  terms  of  the  form 

cos(nt+e,)  ,  exp  [-  Jw  t]  cos  ('^l  .082  -  CZ  u  t  +  e,)  , 

*  S  S  d 

and 

exp[— fu  t]  cos  ("J  0.9241  —  f2  w  t  +  0,)  , 

s  s 

in  which 


(  = 


From  Eq.  (27)  we  would  obtain  terms  of  the  form 

cos(nt  +  O^)  ,  exp(— tWgt)  cosl’J  1  —  J2  w^t  +  0^) 


and 


exp[— fw  t]  w  t  cos ("^  1  —  S2  «  t  +  01) 

So  s  ^ 
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Now,  cos  A  —  cosB  =  2  cos  1/2  (A  —  B);  hence,  the  terms  from  Eq.(26)  could  be  alternatively 
expressed  approximately  as 

cos(nt  +  e^)  ,  exp[— fiUgt]  cos  (^y  1  —  £2  Wgt  +  ep 


and 

exp[-JWgt]  cosNl  -  S2  Wgt  +  0™)  •  sin [W  1.082-  f2-V 0.9241  -  U)  oj^t] 


Now,  sint  for  small  t;  thus,  the  terms  from  Eqs.  (27)  and  (26)  will  be  nearly  the  same  for 
t  <  84  minutes.  This  is  true  regardless  of  the  value  of  the  damping  ratio  f .  Now,  if  C  >  l/2, 
then  by  the  time  that  t  no  longer  approximates  sina)i,eat  quantity  exp[— fwgt]  is  so  small 

that  these  terms  can  be  neglected  anyway.  Thus,  for 

t  <  1/2  and  all  times  t 


or 


we  have 


t  <  84  minutes  and  all  C 


D  «  (S^  +  n^)^  (S^  + 


(28) 


Under  this  assumption  [Eq.  (28)]  the  transfer  functions  of  interest  are  (Fig.  3), 


HAAaA^<S)  =  5 


(S  + 


H 


Zi(t} 


‘A\6V 


(S)  = 


X  R((S  +  rWg)*- +  Wj] 


H 


t..  + 


c2  ^  2 

S  + 


^^*x  CjS(S^  +  R^) 


-“2"8 


7  2  .  2. 

2u)s  <*>2(5  +  w  j  ) 


(S^  +  R^)  (S^  +  Ztu^S  +  ojg)  (S^  +  R^)  (S^  +  ZCUgS  +  Wg  )^ 


t,.  + 


«AA€  <S)  =  2  ‘11 


Cj(S^  +  R^) 


(if.w 

1  s 


-  2  2 
2(v,u;,<a) 
12s 


(S^  +  R^)  (S^  +  2i:w^S  +  uf)  (S^  +  R^)  (S^  +  ZCu  S  + 

8  3  3  8 


t,.  + 


“l"2 


2  2  21 
SC^(S  +  12  )  ^ 


*^1^8 


,  2  2 
Zu.u.u 
12s 


(S^  +  fl^)  (S^  +  2ta>„S  +  (S^  +  R^)  (S^  +  2fu  S  +  oj^)^ 

S  9  s  s 
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H 


AAfiA 


**AXaA^ 


H 


“l“2^22  ^  “1*12 


C^S(S^  +  £2^  S^  + 


2 

Z  8 


Z<i)  Qco^S 
8  Z 


S(S^  +  £2^)  (S^  +  2£a>^S  +  u^)  (S^  +  £2^)  (S^  +  2t«gS  + 


St 


H 


12 


“2*22 


C^(S^  +  £2^) 


in  which 


“2“s 


,  2  _2 
2u  u,S 
8  2 


C^(S^  +  £2^)  (S^  +  2fWgS  +  tOg)  (S^  +  £2^)  (S^  +  2tWgS  + 


H 


“2*12  .  <S^  +  “^  . 


S^  +  £2^  C^S(S^  +  £2^) 


+  <j\) 


2(J-iAs 

8  Z 


C,S(S^  +  £2^)  (S^  +  2t«„s  +  w\)  (S^  +  £2^)  (S^  +  2tu„S  + 

1  8  S  8  8 

2iti. 


"-V®-  T 


XI  '  ‘'1**AA6A„*®' 


s 


s' 


H 


Axav 


(S)  = 


2^^K 


=  C.H, 


«  (S2.2r«  S.a,V  ■ 

s  s 


(S) 


2e“s  = 


2  ,1  ,.2.  2 

“d  =  <*-f  >“s 


We  will  also  need  to  know  the  relationships  among  the  variables  ^  ,  ip  and  ip  and  the  inputs 

X  y  z 

€  ,€  and  €  .  Transforming  Eqs.  (23),  (24)  and  (25)  and  solving,  we  obtain: 

X  y  z 


o2  ^  2 

S  + 


(c  +  :(<  )  + 

'  \r  ^  \rr\' 


“  2  2  '^x  ^xo^  2  2  '^v  "^vo' 

*  S(S  +  £2^)  ^  S  +  n  y 


“l“2 
S(S^  +  12*’) 


(e  +  ^  ) 

'  »  “  wo' 


(29) 


^  = 


,2  ,  „2.  <'x  +  ^xo>  ^  U 


y  (S  +  £2^) 


..  ^~~2  ^*v  *■  '''vo'  ■*  2  ^  2  <^z  ■*  ^zo> 

S  +  £2*^  y  y°  S  +  £2*^  “ 


(30) 


K  = 


S(S  +  £2‘) 


- 


"z  e/oi  .  r.2.  ''x  ^  ■'xo'  ^2  ^  jj2  ''y  ■'yo'  ^  g^g2  ^  jj2j 


=2  2 
S  ■*■  ai  , 

(e  +  ^  )  +  - 5 - 5-  (c  +  ^  ) 

'  V  “vrt'  7  '  Z  ^zo' 


(31) 


in  which  and  ^  denote  the  initial  values  of  ^  ^  and  p  ,  respectively. 

XO  ZO  XV  z 


4 


4! 


I 
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15 


E.  Time-Domain  Responses 

Since  our  ultimate  goal  is  to  find  a  sampled-data  model  to  represent  the  manner  in  which 
the  error  propagates  between  discrete  observation  times,  we  will  need  the  impulse  responses 
associated  with  the  transfer  functions  obtained  in  Sec.  III-D. 

Let 

=  - 2 - 2 - Y - 2 

S(S  +  n  )  (S  +  2ta)„S  +  w“) 
s  s 


^2^®^  ~  2  2  2  ’2  2  ■ 

^  +  il^)  +  2^0)  S+  w  r 

s  s 

We  can  obtain  the  desired  impulse  responses  by  means  of  the  six  transform  pairs  listed 
below: 

Fi(S)  +  I"  cos  nt 

s 


+  exp[- 


'  s 


s'^d 


*.  d 
sinw^t - j  cos 

0) 


“d‘) 


t  >ro 


SF^(S)  *-*•  a  cos  $2t  +  ^  sin  £2t 


/  a-cQoi  V 

+  exp(— JWgt)  (c  cosa)^t+  - -  sinw^tl  t  ^0  , 

'  d  ' 

S^F^(S)  *-►  — aJJ  sin  Ot  +  b  cos  fit 

+  exp[-fa)gt]  |(d- 2ci:Wg)  coswjt-  (cWg(l  -  2t^) 


in  which 


+  cfWg]  sinw^tj 


t  >0 


-2{:[1  +  2(1-  (fi%f)l 


b  = 


1  +  (1-  2f^)  (fi%g) 


4t^  -  1  -  (1  -  lOt^  +  (fiVwg) 
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and 


F^iS)  A  cos  I2t  +  -^  sin  fit  + 


exp[-ta)gtl 

“T^ 

2a)  , 


(®8ina),t+  2  co8a),t 


+  ®  t  sina)jt  +  ®  t  cosa)jt)  t  >  0  , 


SF2(S)  •-*  —An  8in  nt  +  B  cos  nt  + 


exp[-ta)gt) 

- ^ -  (Ssincj^t  +  ©  cosw^t 

2a>  , 


+  ®  t  sina)jt  +  (B)  t  cosojjt)  t  ^  0  , 

2  2  exp[-fa)  t]  ^  ^  ^ 

S  F2(S)  ■"  —An  cos  nt  —  Bn  sin  nt  +  - ^ -  [(— ia)^  ©  +  ®  —  o)^  ®)  sina)^t 

2w  j 

a 

+  (“^^g  ©  +  ©  +  a>^  ©)  cosco^t  +  ®  ©)  t  sinoi^t 


+  (Sofg  ®  +  0)^  ®)  t  cosa)jt]  t  >  0 


in  which 


8a)Jj 


3  +  12S^  -  8f'*  +  (-  5  +  60J^  -  120t 

CO 

s 


+  64fSj 


®=  — y  1+  (3-8S‘)  ^ 

"s  t 

©=  - f  h  +  ^  (2-4f^) 

"s  ^  “s 
®  =  ^  [l  -  2f^  +  ^  (1  -  8f^  +  8f^ 


®  =  J_ 

CO 


-  3  +  2f^  +  (-  15  +  40f^  -  24f^) 


']  ■ 


©  =  - 


[ 


1  +  ^  (2-  12f^) 

CO 

s 


CO 


d 

“Z 

CO 

s 


®  =  -^  I-  1  +  ^  (-  1  +  4f‘) 

CO 

s 


r  O 

=  — ^  1  +  -5 

*^S  0) 


(3-4f^) 
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Using  these  transforms,  we  obtain  the  required  impulse  responses.  Where  possible 
2  2  * 
terms  of  the  order  of  SI  /w  have  been  ignored  with  respect  to  one.  For  certain  values  of  the 

damping  ratio  t,  some  of  the  remaining  terms  in  SI  /to ^  may  also  be  ignored.  The  reader  should 
remember  that,  because  of  the  approximations  made  in  Sec.III-D,  these  expressions,  which  fol¬ 
low,  are  good  for  t  <  84  minutes  and  all  damping  ratios,  or  all  time  and  damping  ratios  greater 
than  l/Z: 


"AXe 


(t) 


— —  cos  S2t  —  sin  fit  +  exp[— tw^t] 


I  —  sin  to  ^t 


to, to 
2  s 


kZ  ,  I  A  1  '3  A  vZ  ,  o  u4, 


Zto^f 


+  ^[(  -  1  +  7r-  lor  +  4f“)  +  cUi-  izr  +  sni - s-  costo  .t 


to, to 
Z  s  , 

+ - r —  t  cos 


2tSl 


“d‘  +  [<^2  -  fl  ^ 


t  >0  , 


^AXc  ^  sin  fit  +  cos  Qt  +  exp(— £a)  t) 
y  s 


s 

g.,  sintojt-costo^t 


tohw 

+  — t  sintojt-  t  cos 

d  tOj 


( 

‘"d‘) 


t  >0  , 


-ZSto, 


V  «i 


^AXe  “  ~  cosnt+  Cj  sin  at  +  exp[-£tOgt]  [  ^3^  |(1  -  f^)  (2J^-  1) 

^d 


+  [(- 1  +  7f  ^  +  4tS  -  (-3  +  12e^  -  8t^] 

"s 


sin  to  . 


7  Z  .  2 

'1  Zt02to.f  to^^i  ? 

+  costo  t  -  -  — -  t  sin  to  t - —  [(l-2f  ) 

d  s 


Zfto^ 

ij 

s 


+  -\  (1  -  8£^  +  8g‘^)J  t  cos 

CO 

s 


‘"d^) 


t>o  , 


^AAe 


f 


2fc|a 


cos  at  +  - —  sin  at  +  exp [—  tto^t] 


\Ci-d 


,  , ,  , ,  sin  to  .t 
s  '  \  u,to  j  d 


—  Q—  costo^t  +  t  sintOjt  — 

1 


j;to  to- 


^j^tsinto^t - costo  t)  t>0  , 

d  1  to^Ci  / 
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-Zfo),  C, 

‘'AAe  “  'cTir'  ^  +  exp[-ta)gt] 

y  1  s  1 


lc,“i  I 


+  (2  -  5t^  -  2f^  +  4t^)|  sinw^t  + 

CJ  ^ 

s 


^‘’  ^ — 77“  COSCJ  .t 

d  d 


ZC^n^S  a),a)g 

- rr- t  sino)  ,t - rt  cosw  .t!  t>0  , 

w_u  ,C.  d  ^  ..  3  d  ‘ 


's^d^i 


Cl“d 


y  y  y  ZO)  y  ^ 

‘^AAe  <t)  “  C2(l  +  2f^  -  n>/)  -  sin  S2t  - 

z  s 


COS  nt 


+  exp[^f<jgt] 


^-f^(5- 4£^)  c,n^(i>  +  4t^) 

- 2 - -  sinWjt—  - 2 -  cosWjjt 


C,(2^  Cy£a^o)  \ 

- t  sincj^t  +  - 2 — ^  t  coswjt  1  t  >  0  , 

“d  o)  ,  / 


“AxaA  <*>  =  mr  exp[-i:a,gtj  sinw^t  t  >0 
X  s 

‘'AXfiV  ^  exp[-tWgt]  sinwjt  t  >  0  . 


‘'AXaA  <-^"s  ®^""d‘  "d^*  coswjt)  t  >  0  , 

y  R«^j 

2 

^Axav  ■  - 3“  “^d^  sinw^jt  +  iTWgW^jt  cosw^t)  t  >0 

y  Ro)  . 


‘'AAaA  ”  C,  ^AxaA^ 
y  lx 


•'AAav  ~  c.  ^Axav 

y  lx 


•^AAaA  ^  C,  ^AxaA 
X  1  y 


‘'AAavJ*^  ■  c.  *^Axav 

X  1  y 


We  will  also  need  the  response  of  the  system  to  initial  conditions.  We  use  the  notation 

^x<‘) 


to  denote  the  response  of  ^^(t)  to  an  initial  position  of  ot  one  unit.  From  Eqs.  (29),  (30)  and 
(31)  of  Sec.III-D,  we  have 


'^x<‘>  2  2 

=  h,  (t)  =  C.S  C“  cos  nt  t  >0  , 

o 
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,(t) 


=  g  (t)  =  Cj,  sin  fit  t  >0 


X  y 


,(t) 


=  h  .  ^  (t)  =  CjC2(l  -  cos  fit)  t  >  0 


r  V  e 

'z  z 

o 


,(t) 


=h^  ,  (t)  =  -C2  8inflt  t>0 


y  X 


Ip  (it) 

—  =  h  .  ^  (t)  =  cos  fit  t  ^  0 


y  y 


>p  (t) 

ir~  ^  ^tP  e  ^  *^1  ^  ° 

r\  •' 


Kit) 

=  h  ,  (t)  =  C,C,(1  -  cos  fit)  t  >  0  , 

o 

>pjn 


=  h  ,  ^  (t)  =  -C^  sin  fit  t  >0 


2  y 


2  2 

f-  =%  e  <‘>  =  <^2  • 

z  ^  z  z 

o 


From  the  fact  that 


^=h^.j,(t)  i,j  =  x,y.z  , 


'Jo  1  J 


and  the  fact  that  ip^,  ip^  and  ip^  appear  only  as  drives  in  the  two  differential  equations  for  AX  and 
AA  [Eqs.(15)  and  (16)],  it  follows  that 

^Mt)  . 

T^=‘'aX€.‘‘’  i  =  x,y.z  , 


and 


AA(t) 


•^AAc.**’  i  = 


Our  time-domain  description  of  the  system  now  lacks  only  the  sixteen  quantities  relating  AX, 
AX,  AA  and  AA  to  ^■^o'  ^^o  ^^o"  transform  Eqs.(15)  and  (16),  including  the 

terms  etc.,  and  solve,  then  by  using  the  approximation  for  D  used  in  Sec.III-D  we  obtain 


=  exp[-{u.gtl  (cosw^t  +  sinwjt)  t  >0  , 

n  n  '  Q  ^ 
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^  =  _L  exp[-Ca,  t]  sin  0)  .t  t  >  0  , 

A\  AA  “d  s  a 

o  o 

u)  C  £2 

expt-fWgt]  (sinw^t- cosw^t)  t  >0 

O  O  OJ^ 


A^  =  exp[-ea)  t]  [-tw  sinw^jt  +  t(a)j  sinw^t 

AA  AX„  £j  , 

o  o  d 

+  tUgdJj  COSUjt)]  t>0 


Also, 

,AX(t).  _  _d  ,AX(t), 
-  dt  '  ■  ’ 


etc. 


F.  Approximate  Time-Domain  Response  for  t  <  84  Minutes 


The  time  responses  given  in  the  preceding  section  (III-E)  are  cumbersome,  to  say  the  least. 
If  we  consider  times  up  to  84  minutes  and  are  willing  to  make  suitable  approximations,  these  time 
responses  can  be  considerably  simplified.  There  are  two  situations  in  which  this  short  time  ap¬ 
proximation  is  applicable: 

(1)  When  the  length  of  time  of  the  whole  operation  does  not  exceed  84  minutes. 

(2)  When  bj;  means  of  external  measurements  each  of  our  seven  variables 

X,  5i,  A,  A,  tpX'  'Py  and  ipz  is  determined  independently  to  an  accuracy  compa¬ 
rable  to  that  which  the  navigator  can  maintain  for  84  minutes;  a  set  of  such 
measurements  is  then  used  at  least  eveiry  84  minutes  to  reset  the  navigator. 

Below  are  listed  the  approximate  time  responses;  those  not  listed  can  be  assumed  to  be  negligi¬ 
ble  with  respect  to  the  given  terms.  These  approximations  have  ignored  terms  of  n/w^{n/w^  » 
l/l7)  or  smaller  with  respect  to  one. 


....  to) 

-jj-Li-  =  1  -  exp[-fWgt]  (cosw^t  +  -^  sinWjjt)  t>0  , 

y  Vo  ^ 


•^AXe  (*)  = 


‘’AXfiV  t  >0  , 


^AX6A  <*>  =  mr  exp[-e«gt]  sina,jt  t  >0 
X  s 


‘^AAc,<t>  = 


^AA6V  -  ^  exp[-Sa,gtJ  sinw^t  t  >0 

•'  1 


‘^AAdA  =  -R^  sin"d‘  ^ 

y  is 


\  e  = 

X 


,(t) 


1  t  >  0  , 
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Ip  (t) 

=  =  l  t>0  . 

y  y 


,(t) 


%  c  <*)  =  . 
z  z^ 


=  1  t  ^0  , 


=  exp[-fa)gt]  (costJ^t  +  -jj-2.  sinWjjt)  t>0 
*0  d 


—  =  i;r~  exp[-£ajgt]  Sinoj^t  t  >0 

AX  “d 
o 

AA(t)  AX(t)  AA(t)  _  AX(t) 

AA  “  AX  •, 

o  o  AA  AX 

o  o 


AX(t)  _  AA(t)  d  .AX(t).  AX(t)  AA(t)  _  d  .AX(t), 

A>  AA  ~  Ht  Iax  ‘  ■  v  ■  dt  *  •  ‘ 


AX  AA  '  dt  ^AX 
o  o  o 


AX  AA 
o  o 


AX. 


G.  Equivalent  Discrete-Time  System 


Using  the  results  of  Sec.  III-E  or  III-F  (whichever  set  of  responses  is  applicable),  we  can 
now  write  the  set  of  matrix  equations  which  describes  the  error  propagation  between  sample  (ob¬ 
servation)  times.  Let  x(t)  denote  the  seven-tuple  column  vector: 


Then, 


in  which 


AX(t)  =  x^ 
AX(t)  =  Xj 
AA(t)  =  Xj 
AA(t)  =  x^ 

yt)  =  Xg 
^^y(t)  =  X^ 


7 

x(t  +  T)  =  ♦.^.^xlt)  +  2  xNt  +  T)  , 
i=l 


(32) 
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A\(T)  AX(T) 

^^0  AX 

o 

AX(T) 

d  AXftl 
dt  AX 

o 

etc. 

t=T 

AA(T)  AA(T) 

'■  AX  ... 

0  AX^ 

AA(T) 

-KK- 

o 

AA(T) 

'•aa,  ITI 
y 

‘‘AAe.^T) 

*77  = 

d  AA(t) 
dt  AX 

0 

etc. 

t=T 

0 

0 

0 

0 

%  e  <T) 

’^x  y 

^  c  <T) 

^X  z 

0 

0 

0 

0 

N  , 

X 

\  c  <T) 
y  y 

y  z 

0 

0 

0 

0 

\  c  <T) 

X 

\  e  <T) 
z  y 

%  e  <T) 
z 

and  the  term  which  appears  in  the  k*^  row 

of  +  T) 

is 

xj 

*k 

T 

=  r  h  (T)  y  (t  +  T  - 

Jo  Vj  J 

-  T)  dr  , 

j,k=l,2 . 7  , 

yi  =  =  6A^.  yj  =  6Vy,  y^  =  6Ay,  yg  =  €^,  y^  =  y^  =  . 

It  should  be  noted  that  the  second  and  fourth  rows  of  are  the  derivatives  of  the  first  and  third 
rows,  respectively,  the  derivatives  being  evaluated  at  t  =  T  as  indicated.  Similarly,  we  note  that 

. '  • 


Any  impulses  appearing  in  these  derivatives  must  be  retained  in  the  impulse  responses  h^-^^y  (t) 
and  h.'.  (t) .  ^ 

Our  ultimate  goal  is  to  be  able  to  study  the  error  propagation  when  we  couple  our  inertial 
navigator  to  external  measurements  in  some  optimum  fashion.  This  problem  seems  to  admit  of 
treatment  most  readily  when  viewed  in  terms  of  Kalman's  formulation  of  the  optimum  filter  prob¬ 
lem.  In  order  to  make  our  model  amenable  to  treatment  via  Kalman's  method,  we  must  be  able 
to  express  the  inputs  x^(t  +  T)  as  the  responses  of  linear  (discrete)  systems  to  white  noise.  Since 
our  estimation  procedure  will  be  based  strictly  on  first-order  statistics,  we  must  ensure  only  that 
the  matrices 


R.  .(NT)  =  N  =  0,  ±1,  ±2,  .  . . 

i.j  =  1.2 . 7 
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! 


are  the  same  for  both  our  original  inputs  and  the  responses  to  the  linear  systems.  We  will  assume 

that  the  inputs  e  (t) ,  e  (t),  £_(t),  6V  (t),  6A  (t),  dV  (t)  and  dA  (t)  are  all  mutually  statistically 
X  y  z  X  X  y  y 

independent.  Thus  we  need  consider  the  matrices  R..  only  for  i  =  j,  i  =  1,  2,  ....  7. 

5  6  7  ” 

Consider  first  x  ,  x  ,  x  .  We  assume  that 


g  (t)  =  E{e.(t)e.(t  +  T)}  =  expf-v^lTl]  i  =  x,y, ! 


,.th 


Thus,  the  j-k  entry  of  R  ^  (NT)  is 

^i  i 

=  I  I  V  £.‘^2)  ♦£.£.<NT  +  <Pi  -  <Pi} 

”  *'0  11  K  i  II 


''£.  »X.£.Hx,£. 
1  3  1  k  1 


N  =  0 


in  which 


“x  £  “x  £  exp(- vJntI  )  N  >  0 
i  31  k  i 

<"£  «x£  “x  £  ^^Pl-^ilNTll  N<0  , 

£i  x^£^  x^£.  1 


nT 

^  h^^£_(,>2)  exp  [-^il  ^1  -  ,,2!  1 

+  r'^ 

«x.£.=J^  expt+,/jV>ildv^  , 

«x.£.=3^  exp(-...V,jld^^  , 


i  =  x,y,  z  3,  k  =  1,2, ..  .,7 
x^  =  AX,  x^  =  AX,  Xj  =  AA,  =  AA,  Xj  =  x^  = 

Now,  consider  the  representation 

<t)  +  B..d,(t-T)l  , 


3« 


31  1' 


in  which 


dj(t)  =  ajd.(t  -  T)  +  w.(t) 


di(t)  =  21  Wj(t  -  mT)  (a.)*”  i=x,y,  z 

m=0 


(33) 


(34) 


(35) 


(36) 


(37) 


(38) 


(39) 


The  w's  are  all  white  random  processes  with  zero  mean,  and  the  w^'s  are  taken  to  have  unit  var¬ 
iance.  The  variance  of  the  w,  's  and  the  correlation  between  all  the  w's  are  left  free  for  the  time 

J*i 

being. 
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Substituting  Eq.  (39)  into  Eq.  (38)  and  noting  that 


Wj(t)Wj(t  +  NT)  =  WjWja(NT) 

we  obtain 


(40) 


"i . "i, 


X.  *(t)x^‘(t  +  NT)  =  r.yNT)  , 

(w.  w, —  +  - .)  N  =  0 

,J"i^‘i  l-a.2  31  ki/ 


N  >  0 


N  <  0 


Our  discrete  white-noise-generated  inputs  will  then  match  our  original  inputs  if 


a.  =e 


-i/.T 
1 


H.  H,  =  w.  w,  +  - i— 7  B..B,  . 

je.  ke.|  j€.  kc.  jx  ki 


B..B,  . 

a r  H7,  H  ■  =  B.,ar  *  w,  ,  w.  +  — J---  -- 


.  -1: 


*Xi  —  vf.x  vr .  T  «  . 

i  3"i  ‘‘"i  3^  ^  “"i  ‘ 


(41) 


(42) 

(43) 

(44) 


(45) 


The  solutions  to  Eqs.  (43),  (44)  and  (45)  are 


w.w.  =  IH.  — 


1  J€.  €.  \  J€.  .  Z 

J  1  i  \  ^  I  1  —  a. 


H.  H. - 

]€.  ke.  .  „2  I 

■'  1  i  1  —  a.  / 


i  =  X,  y,  z  j,  k  =  1,  2,  .  .  .  ,7 


(46) 


(47) 


(48) 


as  can  easily  be  verified  by  substitution.  We  can  thus  obtain  equivalent  discrete  inputs  for  the 
signals  €^(t),  €y(t)  and  e^(t). 

We  now  consider  the  inputs  6V  (t),  6A  (t),  6V  (t)  and  6A  (t).  Our  analysis  here  will  be 

X  X  y  y 

much  easier,  as  we  can  assume  that 
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U-gg’4MH 


Fig.  4.  Summary  of  diicr«to-time  error  propogation.  The  comporrenfs  of  4^  and  fhe  term* 

”ik'  ^fk'  l^im^kml 


2i 


<l>  (t) 

yi^i 


(p.  exp[-t/.|TlJ 


i  =  1,  2.  3,4 


6V^,  72  =  «A^.  73  =  6V  .  74  =  «A 


where  is  large  enough  so  that 


1 


V 


(t)  7j(t  -  t)  dr 


J 


is  uncorrelated  with 


1 


h„  „  (t)  7.(t  +  NT  -  T)  dr 

’‘kyj  J 


for 


N  ^  0  ,  j  =  1, 2,3.4,  k=  1,2.. ...7  . 

Thus,  these  four  inputs  can  be  replaced  b7  the  discrete  white  inputs 


kj 


1,2,  ....7;  j  =  1,2.3.4 
^T  rT 


'^kj^nm  =  \  So  So  «P  I-'kl  1 1  '^’’l'^<^2 


where 


where 


m  =  j  n  =  k 

=  0  otherwise 

k  =  1,  2 _ ,7 

j  =  1.2,  3,  4 

Recall  that  we  have  assumed  that  6V  ,  6V  ,  6A  and  6A  were  statisticall7  independent.  In  prac- 

X  y  X  y 

tice  fiV^  and  dV^  will  not  be  uncorrelated.  This  can  easil7  be  corrected  b7  setting 


"  S  S^  Vfiv  <’’!>  •'x.dv  <^2>  Rfiv  6v  <'^1  -  ’’2)  ‘‘’’l‘^<^2 

If  V  IV  V  V 


k  X  "j'  7  ■  X'  y 

j,  k  =  1,  2,  .  .  . ,  7 

We  now  have  all  the  information  needed  for  our  discrete-time  model.  This  information  is 

summarized  in  Fig.  4.  For  operating  times  less  than  84  minutes  the  impulse  responses  of 

Sec.III-F  ma7  be  used;  otherwise,  those  of  Sec.III-E  must  be  used.  For  picket  ship  operation 

the  components  of  and  the  terms  H.^  ,  H.”  and  |H.  H,  _  I  are  constant.  For  normal  ma- 
^  77  jm  ]m  '  jm  km' 

rine  operation  the  coefficients  of  the  differential  equations  governing  error  propagation  (15, 16, 

23,  24  and  25)  var7  so  little  during  a  four-hour  period  that  the7  can  often  be  considered  constant. 

Hence,  for  T  <  4  hours,  and  the  quantities  H.t_,  H.”  and  1h._H, _ I  can  often  be  determined 

77  ’  jm  jm  '  jm  km 

at  each  sample  time  b7  using  the  given  expressions  with  the  appropriate  values  of  and  in¬ 
serted  into  the  impulse  responses. 
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IV.  ESTIMATION  OF  NAVIGATOR  ERROR 


A.  Estimation  Equations  When  Navigator  Is  Not  Reset  Following  Each  Observation 

In  developing  our  solution  to  the  estimation  problem  we  will  make  repeated  use  of  two  prop¬ 
erties: 

(1)  Let  X|  and  X2  be  two  random  variables,  and  let  2^  and  X2  be  the  best  lin¬ 
ear  estimates  of  x^  and  X2,  respectively,  using  the  same  observed  data 
for  each  estimate.  Then,  the  best  linear  estimate  of  ax.  +  bx2,  again 
using  the  same  data,  is 

ax^T^2  =  3.2^  +  b22 

(2)  Consider  any  fixed  class  F  of  functionals  on  some  given  data  with  the 
property  that  any  linear  combination  of  functionals  in  F  is  again  in  F. 

Let  X  be  the  functional  in  this  class  that  is  the  best  estimator  of  x. 

Then,  (x  —  2)  is  uncorrelated  with  any  functional  of  the  given  class, 
that  is, 

E{(x-  2)  f}  =  0  fcF,  2€F 

In  the  above  statements  "best"  is  used  in  the  minimum-mean- square  error  sense.  Property  (2) 
follows  directly  if  we  regard  mean-square  estimation  as  finding  a  projection  in  an  inner  product 
space.  However,  for  completeness,  simple  proofs  of  both  properties  are  given  in  the  appendix. 

For  convenience  we  take  T,  the  time  between  observations,  to  be  one  unit.  Then  the  model 
for  our  process  is 

z(t)  =  H(t)x(t)  +  v(t)  ,  (49) 

x<t)  =  ♦(t)x(t-  1)  +  w(t-  1)  ,  (50) 

in  which  v(t)  and  w(t)  are  white  and  have  zero  mean.  Since  it  suits  our  present  purpose,  we  also 
assume  that  w  and  v  are  uncorrelated.  Now  our  estimator  2^  is  linear  and  may  be  written 

00 

2(t)  =  jA(k)z(t-k)  .  (51) 

k=0 

The  sum  in  Eq.  (51)  actually  extends  back  only  over  all  past  values  of  data;  we  use  *>  as  the  up¬ 
per  index  for  convenience.  We  now  define 

00 

f 

z*(t)  =  z(t)-  Yj  B^z(t-k) 
k=l 

rr* 

such  that  E  =  0;  k  =  1, 2, . . .  .  We  can  then  express  2(t)  as 

oe 

2(t)  =  C(t)  af(t)  +  Y  •  (52) 

k=l 

Now,  because  z’*'(t)  is  uncorrelated  with  all  the  z(t  —  k),  k  =  1,  2,  .  .  .  ,  the  minimization  equations 
separate  into  two  sets  of  equations,  one  involving  only  C  and  the  other  involving  only  the  C^^'s. 
Thus,  the  Cj^'s  can  be  determined  independently  of  C  and 
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Yj  —  k)  =  best  linear  estimate  of  x^., 

given  z(t—k),  k  =  1,  2,  .  .  .  ; 

=  best  linear  estimate  of  *(t)3c^_^  + 
given  z(t  —  k),  k  =  1,  2, .  .  .  ; 

=  *(t)fi(t  -  1)  +  0  . 

Equation  (53)  follows  directly  from  property  (1).  Thus, 
x(t)  =  C(t)z*(t)  +  ♦(t)£(t- 1) 

Now,  z*(t)  is  uncorrelated  with  *(t)2(t—  1); 

Z*(t)  =  z(t)-B(t)*(t)S(t- 1)  : 

and 

0  = 

=  2^^J^*'^(t)  -  B(t)*(t)  x^_^^'^l*(t)’^ 

We  assume  that  i{t)  and  are  invertible;  thus, 

B(t)  =  it- A-1  ^  • 

Now,  from  Eqs.(49)  and  (50) 

2t^-l  ■  ^-t-1  -t-t-1  "  -t^-1 

=  H(t)  [♦(t) 

=  H(t)  ♦(t)  . 

and  from  property  (2)  we  have 

or 

T  O  ^  T 

Substituting  Eqs.(57)  and  (56)  into  Eq.  (55),  we  obtain 
B(t)  =  H(t)  ♦(t)  ^t-A-1  ^ 

and 

»(t)  =  C(t)  (z(t)  -  H(t)  ♦(t)  ft(t-  1)]  +  *(t)  l(t  -  1) 


(53) 


(54) 


(55) 


(56) 


(57) 


(58) 


(59) 
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We  now  need  to  solve  for  C{t) .  From  property  (2)  we  have 


=0 

or 


{Xt  -  C(t)[Zt  -  H(t)  ♦(t)  -  ♦(t)  x^_i  }  zT"  =  0  ' 

or 

-  *(t)  =  C(t)[zjZ^'^  -  H<t)  ♦(t)  • 

Now,  by  Eqs.(49)  and  (50),  we  have 

A  T  A  T  „T,.,  ,  7  T 

Xt.lZt  =  Xt.iXt  H  (t)  +  Xt.iVj 


A  T  itT,., 


=  Xj.iXtTi 

We  can  express  the  covariance  matrix  of  the  error  as 
S(t  -  1)  =  (x^.i  -  k^.^)  (XtTi  -  StTi  > 


=  (Xj.!  -  Xt^i  -  -  2j,i) 


and  thus 


^t-l^t-l  =  2it,iXtIi  -S(t-l)^X(t-l)-S(t-l) 

and 


=  X(t  -  1)  ♦'^(t)  H''^(t)  -  S(t-  t)  ♦'^(t)  H'^(t) 
From  Eq.  (49)  we  also  obtain 


(60) 


(61) 


(62) 


(63) 


Also, 


=  H(t)  x^'^H'^(t)  +  +  H(t) 

=  H(t)  X(t)  H'^(t)  +  V(t)  ,  (65) 

T 

in  which  V{t)  denotes  .  Further, 

X(t)  =  =  *(t)  +  *{t)  Xt.jWtTi 

=  *(1)  X(t  -  1)  *'^(t)  +  W(t  -  1)  ,  (66) 

in  which, W(t)  denotes  •  Substituting  Eqs.  (63)  through  (66)  into  Eq.  (60),  we  obtain 

♦(t)  X(t  -  i)  *'^(t)  H’^(t)  +  W(t  -  1)  H'^(t)  -  *(t)  X(t  -  1)  ♦'^(t)  H’’''(t)  +  ♦(t) 

X  L(t-  1)  ♦'^(t)  H'^(t)  =  C(t)  [H(t)  *(1)  X(t-  1)  ♦’^(t)  H'^(t)  +  H(t)  W(t-  1) 

X  H'^(t)  +  V(t)  -  H(t)  ♦(t)  X(t  -  1)  ♦'^(t)  H'^(t)  +  H(t)  ♦(t)  S(t  -  1)  ♦'^(t)  H’^(t)] 

and 

C(t)  =  (W(t-  1)  +  *(t)  2(t-  1)  ♦'^(t)]  H'^(t){V(t)  +  H(t)  [W(t-  1) 

+  ♦(t)  S(t  -  1)  ♦'^(t)]  H'^(t)}‘*  (67) 

Note  that  C(t)  depends  on  r(t  -  1);  thus,  in  order  to  complete  our  recursive  scheme,  we  must  be 
able  to  express  Z(t)  in  terms  of  C(t)  and  i:(t—  1): 

2(t)  =  (Xj  -  x^)  {x^  -  ^  I  =  -  S^^) 

=  X(t)  -  x^Zt"^  c'^(t)  +  C^(t)  -  ♦'^(t)  (68) 

Using  Eqs.  (64),  (56)  and  (57),  we  obtain 

Z(t)  =  X(t)  [I  -  H'^(t)  c'^(t)l  -  ♦(t)£j_jX^'^j  ♦'^(t)  (I  -  H’^(t)  c’^(t))  ,  (69) 

and  using  (62)  and  (66),  we  obtain 

E(t)  =  {  *{t)  X(t  -  1)  ♦(t)’^  +  W(t  -  1)  -  ♦(t)  [X(t  -  1) 

-  r(t  -  1)1  *'^(t)}  [I  -  H'^(t)  c'^(t)]  , 

S(t)  =  [W(t  -  1)  +  ♦(t)  2(t  -  1)  ♦'^(t)]  (I  -  H'^(t)  c'^(t)]  (70) 

B.  Estimation  Equations  When  Navigator  Is  Reset  Following  Each  Observation 

The  equations  in  Sec.IV-A  were  developed  under  the  assumption  that  the  estimate  of  the  nav¬ 
igator  error  was  not  used  to  reset  the  quantities  in  the  navigator.  Here  we  assume  that,  follow¬ 
ing  an  observation,  x(t)  is  used  to  reset  the  navigator  (either  by  resetting  the  navigator  computer 
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or  both  the  computer  and  stable  platform) .  If  we  use  ^(t  —  1)  to  reset  the  navigator,  then  the 
model  for  our  error  propagation  becomes 

z'(t)  =  H(t)x'(t)  +  Y(t)  ,  (71) 

x'(t)  =  *(t)[  x'(t- l)-i'(t- 1)1  +  w(t- 1)  ,  (72) 

in  which  we  have  used  primes  to  distinguish  the  quantities  in  reset  operation  from  those  in  nor¬ 
mal  operation.  The  relations  between  these  two  sets  of  quantities  are; 


x'(0)  =  x(0)  x'(0)  =  x(0)  =  0 

X'(l)=x(l)  x'(l)=x(l) 

x'(2)  =  ♦2[x(l)  -  ^(1)1  +  w(l)  =  x(2)  -  *(2)  gii) 

Thus, 

i<{2)  =  $(2)  -  *(2)i(l) 

and 

x'(3)  =  *(3)  [  x'(2)  -  x'(2)  +  w(2) 

=  *(3)  [x(2) -$(2)1  +  w(2)  =x(3) -*(3)8(2) 

Thus, 

8'(3)  =£(3)-*(3)8(2)  . 

In  general. 

X'(n)  =  x(n)  —  ♦(n)x(n— 1)  , 

2'(n)  =  x(n)  —  ♦(n)x(n—  1) 

and 

z'(n)  =  H(n)  x'(n)  +  V(n)  =  z(n)  -  H(n)  *{n)  8(n  -  1)  . 

Substituting  into 

x(n)  =  C{n)  [  z(n)  —  H(n)  ♦(n)£(n—  1))  +  *(0)8(0  -  1)  , 

we  obtain 

2'(n)  +  *(n)8(n-l)  =  C(n)[z'(n)  +  H(n)  *(n)  8(n  -  1)  -  H(n)  *(n)  x(n 

+  *(0)  R(n  -  1) 


or 


5'(n)  =  C(n)  2i(n) 

The  quantity  C(n)  is  obtained  exactly  as  before,  since 


z:-(n)  =  (x;-g;)(x;-g;)'^ 


(73) 


=  <^n-^n>  <^n-^>^  ’ 

The  estimation  scheme  thus  remains  basically  unchanged.  The  error  is  unchanged,  as  is  the 
basic  computational  complexity,  in  that  we  must  still  generate  the  matrices  C(t)  and  £(t)  in  the 
same  recursive  fashion.  The  estimator  Eq.  (73)  is  slightly  simplified  over  Eq.  (58). 
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C.  Summary  and  General  Remarks 

In  summary,  we  have  the  two  possible  models 


z(t)  =  H(t)  x(t)  +  y(t) 

(Eq.(49)] 

x(t)  =  *(t)x(t- 1)  +  w(t- 1)  , 

(Eq.(50)] 

or 

z'(t)  =  H(t)x'(t)  +  y(t) 

[Eq.(71)] 

S'(t)  =  *(t)  [x'(t-  1)  -S'(t-  1)1  +  W(t-  1)  . 

[Eq.(72)] 

The  two  estimators  are  given  by 

i(t)  =  C(t)  [  z(t)  -  H(t)  *{t)  g(t  -  1)1  +  *(t)  g{t  -  1) 

[Eq.(59)] 

and 

x'(t)  =  C(t)  z(t)  . 

[Eq.(73)] 

The  matrix  C(t)  is  determined  by  the  pair  of  recursive  relations 

C(t)  =  [W(t  -  1)  +  ♦(t)  S(t  -  1)  *'^(t)]  H’^(t)  {V{t)  +  H(t)  [W(t  -  1) 


+  ♦(t)  S(t-  1)  ♦'^(t)l  H'^(t)}'^  [Eq'.(67)] 

S(t)  =  [W(t-  1)  +  ♦(t)S(t-l)*'^{t)]  [I-H'^(t)c'^(t)]  .  [Eq.(70)] 

in  which 

2(t)  =  covariance  matrix  of  the  estimation  error  x^  — 

W(t  —  1)  =  covariance  matrix  of  the  input  w(t  —  1) 

V(t)  =  covariance  matrix  of  the  measurement  error  v(t) 

Our  computation  starts  with  some  value  for  S(o).  If  the  platform  were  initially  perfectly 
oriented  and  the  initial  position  known  exactly,  then  S(o)  would  be  zero.  This  will  not  usually 
be  the  case,  and  S(o)  will  have  to;  be  determined  from  the  alignment  procedure  used.  Knowing 
S(o),  we  calculate  C(l)  from  S(o)  via  Eq.  (67).  Then  2(1)  can  be  calculated  from  C(l)  via  Eq.  (70). 
This  process  can  be  repeated  continuously  and  we  obtain  2(t)  for  t  =  1,  2,  .  .  .  ;  in  the  process  we 
have  determined  C(t)  and  hence  $(t).  If  ♦(!)  is  constant,  i.e.,  ♦(t)  5  ♦,  then  the  steady- state  value 
of  2  can  possibly  be  obtained  by  setting  2(t)  =  2(t  —  1)  and  eliminating  C(t)  between  Eqs.(67)  and 
(70). 

This  discussion  concerning  the  steady-state  value  of  2(t)  naturally  poses  the  question;  Under 
what  conditions  is  the  steady- state  value  of  2(t)  unique?  Another  point  that  requires  an  answer 
is:  When  is  the  optimum  estimator  asymptotically  stable?  The  asymptotic  stability  of  the  esti¬ 
mator  assumes  importance  because  we  do  not  wish  small  bias  inputs,  which  have  been  neglected 
in  our  analysis,  to  be  able  to  cause  arbitrarily  large  errors  in  our  steady-state  estimate.  To 
these  questions  Kalman®  provides  the  following  answer,  which  is  definitive  but  perhaps  somewhat 
cumbersome  to  apply. 

Theorem  (Kalman): 

Let  the  system  defined  by  Eqs.  (1)  and  (2)  be  completely  observable  and  completely  control¬ 
lable  .  Then 
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(1)  The  optimum  estimator  is  uniformly  asymptotically  stable. 

(2)  All  solutions  corresponding  to  different  choices  of  a  covariance  matrix  S(o) 
converge  uniformly  to  the  same  solution. 

Kalman  refers  to  the  system  defined  by  Eqs.(l)  and  (2)  as  completely  observable  if  for  every 

T 

vector  P  and  every  time  t^  there  exists  a  Tit^)  and  an  unbiased  estimator  n  of  P  x(T),  in  which 
fi  is  a  linear  function  of  z(t),  tjj  •$  t  <  T.  He  calls  this  same  system  completely  controllable  if 
there  exists  a  forcing  function  w(t)  which  can  take  the  system  from  rest  to  any  desired  state  in  a 
finite  time.  For  mathematically  definitive  conditions  for  complete  observability  and  complete 
controllability  and  a  thorough  discussion  of  the  whole  question,  the  reader  is  referred  to  Chap¬ 
ters  15  and  16  of  Ref.  5. 

One  further  point  deserves  mention.  The  estimator  described  gives  simultaneously  the 
minimum-mean- square  error  estimate  of  each  component  of  x.  Our  concern  may  center  on  the 
position  estimate  in  some  particular  direction,  say  in  the  quantity  x^  cos  6  +  x^  sin  6,  and  we 
may  wonder  if  it  is  possible  to  find  another  estimator  which  gives  a  better  estimate  of  cos  6  x^  + 
sin  0  However,  the  estimate  obtained  using  our  estimator  is  cos  0  +  cos  0  which,  ac¬ 

cording  to  property  (1),  is  the  minimum-mean- square  estimate  of  cos  0  x^  +  sin0  x^  for  the  given 
data.  This  property  of  our  estimator  is  equivalent  to  a  statement  that  the  ellipsoid  of  concentra¬ 
tion  of  our  estimation  error  lies  wholly  within  the  ellipsoid  of  concentration  obtainable  by  any 
other  linear  estimator.^ 

V.  COMPUTATIONAL  ASPECTS 

In  order  to  use  the  formulation  presented  in  Sec.  IV  to  find  the  optimum  estimator  or  to  solve 
for  the  minimum  possible  estimation  error,  considerable  machine  computation  is  necessary  in 
almost  all  cases  of  interest.  Some  of  the  computations  can  be  carried  out  on  an  analog  computer, 
but  the  use  of  a  digital  machine  at  some  points  is  unavoidable.  We  review  here  the  basic  quanti¬ 
ties  that  must  be  calculated. 

In  order  to  apply  the  formulas  given  in  Sec.  IV  there  are  five  basic  matrices  which  must  be 
calculated:  ♦(!),  W(t),  H(t),  V(t)  and  2(o).  The  first  two  of  these  are  found  from  the  description 
of  the  inertial  navigator,  the  second  two  from  the  description  of  the  measurement  process  and  the 
last  from  the  description  of  the  initial  alignment  process. 

In  obtaining  the  matrices  ♦(t)  and  W(t)  there  are  98  quantities  which  must  be  found.  Forty- 
nine  of  these  are  the  response  at  time  T  of  each  of  the  seven  basic  error  quantities  (AA,  AA, 

AA,  AA,  >l>  ,  Ip  ,  ^  )  to  an  initial  unit  displacement  of  each  one  of  these  same  seven  quantities. 

X  y  z 

The  second  49  quantities  required  are  the  responses  at  time  T  of  each  of  the  above  seven  quan¬ 
tities  to  a  "unit  impulse"  excitation  of  each  of  the  seven  inputs  (6V  ,  6A  ,  6V  ,  6A  ,  e  ,  €  ,  c  ) . 

X  X  y  y  X  y  z 

These  98  quantities  are  given  in  Sec.III-E  for  all  T.  An  approximate  set,  which  can  be  used  when 
the  time  of  the  total  operation  does  not  exceed  84  minutes,  is  given  in  Sec.III-F.  Fortunately,  not 
all  of  these  98  quantities  are  distinct;  as  pointed  out  in  Sec.III-E,  many  of  the  quantities  in  the  set 
of  49  initial  condition  responses  are  identical  to  quantities  in  the  set  of  49  impulse  responses.  If, 
for  some  reason,  it  is  not  desirable  to  use  the  responses  given  in  Sec.III-E,  they  can  be  calcu¬ 
lated  in  two  ways.  The  first  way  would  be  to  use  a  digital  computer  to  numerically  integrate  out 
the  responses  in  accordance  with  Eqs.(15),  (16),  (23),  (24)  and  (25).  An  easier  way  might  be  to 
use  an  analog  computer  to  simulate  these  five  equations  and  find  the  desired  responses  by  making 
the  necessary  number  of  runs  on  the  analog  computer.  In  finding  the  impulse  responses,  any  in¬ 
put  of  duration  less  than  l/20^^  of  the  smallest  time  constant  of  the  system  would  be  suitable. 
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As  shown  in  Fig.  4,  *(t)  consists  of  a  7  x  7  matrix  plus  several  smaller  arrays.  The 

expression  for  4^^  is  given  in  Sec.III-E  in  terms  of  the  initial  condition  responses  mentioned 

above.  Three  of  the  remaining  terms  in  ♦(!)  are  the  time  constants  of  the  correlation  functions 

+  1 

of  the  gyro  drifts  and  must  be  estimated.  The  other  quantities  in  *(t)  are  the  H  s.  These  must 
be  calculated  by  the  integration  in  Fig.  4.  This  integration  could  be  carried  out  by  means  of  ei¬ 
ther  digital  or  analog  equipment;  it  can  also  be  performed  analytically  by  hand,  with  perhaps  a 
reasonable  amount  of  patience. 

As  shown  in  Fig.  4,  the  entries  of  the  W(t)  matrix  are  found  by  calculating  the  quantities  H^, 
H"  and  .  The  quantities  H”,  like  the  quantities  are  found  by  a  single  integration 

involving  an  impulse  response,  as  indicated  in  Fig.  4.  The  quantities  involve  a  double 

integration  involving  two  impulse  responses  simultaneously.  This  double  integration  can  be  car¬ 
ried  out  on  a  digital  machine  and,  perhaps  somewhat  awkwardly,  on  analog  equipment.  It  can 
also  be  carried  out  by  hand,  but  the  task  is  so  unwieldy  that  it  would  require  an  extreme  amount 
of  patience  and  devotion  to  duty. 

We  can  only  indicate  here  how  the  matrix  H(t)  is  obtained,  since  we  have  not  specified  the 
method  of  taking  measurements.  We  assumed  that  the  measurements  could  be  expressed  as 


r(t)  =  F[q(t)  +  x(t)]  +  v(t)  ,  (74) 

in  which  q(t)  represents  the  true  value  of  the  seven  quantities  used  to  describe  the  navigator  op¬ 
eration.  The  function  F  must  be  found  from  the  geometry  of  the  situation.  The  matrix  H(t)  is 
then  expressible  as 


Hi.(t) 


aF.(s(t)) 

8s.(t) 


s(t)  =  q(t) 


(75) 


Equation  (75)  can  be  evaluated  analytically  or  by  a  numerical  differencing  approximation.  The 
matrix  V(t)  must  be  found  by  an  enlightened  examination  of  the  errors  involved  in  the  measure¬ 
ments. 

The  matrix  2(o)  can  be  estimated  only  by  someone  well  versed  in  the  procedures  used  to 
align  an  inertial  platform.  For  an  introduction  into  the’ many  possible  means  of  alignment,  the 
reader  is  referred  to  Chapter  6  of  Ref.  4. 

With  the  five  matrices  ♦(t),  W(t),  H(t),  V(t)  and  S(o)  available,  the  recursive  formulas  sum¬ 
marized  in  Sec.IV-C  can  be  employed.  However,  it  should  be  mentioned  that  these  formulas  can 
be  written  in  various  ways  and  the  best  choice  depends  on  the  problem.  The  following  is  a  rear¬ 
rangement  of  the  basic  formula  that  is  not  self-evident.  Define 


I(t)  =  S'^t)  ;  (76) 

then,  in  place  of  Eqs.(67)  and  (70)  we  have 

I(t)  =  [  *(1)  I’^t  -  1)  ♦’^(t)  +  W(t  -  1)  ]  "S  H’^(t)  V"Vt)  H(t)  (77) 

For  the  actual  estimator,  Eq.  (59)  can  be  replaced  by 

S(t)  =  I"^t){[*(t)r^t- l)*'^(t)  +  W(t- !))■*  ♦(t)»(t-l) 

+  H'^(t)V^t)z(t)}  .  (78) 
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The  necessary  matrix  manipulations  needed  to  obtain  Eqs.  (77)  and  (78)  are  outlined  in  Statement  3, 
Appendix  A.t 

The  formulas  of  Sec.  IV  require,  at  each  step,  the  inversion  of  a  matrix  whose  rank  equals 
the  number  of  measurements  made  (the  number  of  components  of  z).  Equations  (77)  and  (78)  on 
the  other  hand,  require  the  inversion  of  a  matrix  with  rank  equal  to  the  number  of  states  (the 
number  of  components  of  x) .  In  either  case,  matrix  inversion  is  required  and  there  are  many 
possible  techniques  available.  Reference  7  is  one  of  many  sources. 

In  the  case  in  which  t,  W,  H  and  V  are  constant,  the  steady-state  value  of  Z  is  the  solution 
of  the  equation 

S  =  F  -  Fh'^[V  4  HFh'^]"*  HF  ,  (79) 

in  which 

F  =  W  4- 

This  matrix  equation  results  in  a  system  of  simultaneous  quadratic  equations  which  may  be  solv¬ 
able  by  machine  methods. 

There  is  one  potentially  important  technique  for  simplifying  each  recursive  computation  con¬ 
siderably  while  possibly  increasing  the  number  of  recursive  calculations.  Recall  that  the  damp¬ 
ing  ratio  f  is  zero  if  the  velocity  log  is  not  incorporated  into  the  continuous  system.  With  f  =  0, 
the  expressions  for  H  ,  H”,  entries  of  simplify  by  almost  an  order  of  mag¬ 

nitude.  The  information  given  by  the  velocity  log  can  still  be  used  in  correcting  position  errors 
by  sampling  the  velocity-log  signals  at  a  rate  twice  the  "bandwidth"  of  these  signals.  If  the  sam¬ 
pling  time  T  of  our  system  is  set  equal  to  the  time  between  velocity- log  samples,  then  these  sig¬ 
nals  can  be  incorporated  in  2(t)  and  used  in  an  optimum  manner  to  help  obtain  the  estimate  of  j((t). 
This  procedure  may  also  increase  the  over-all  accuracy  of  the  system,  since  the  conventional 
velocity-log  damping  does  not  necessarily  employ  the  velocity  data  in  an  optimum  manner. 


tThe  formulas  con  olio  be  derived  directly  by  using  o  line  of  development  different  from  thot  in  Sec.  IV. 
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APPENDIX  A 

DERIVATION  OF  CERTAIN  PROPERTIES  OF  THE  OPTIMUM  ESTIMATE 
AND  AN  ALTERNATE  FORM  FOR  THE  ESTIMATOR 

Statement  1:—  Consider  the  data  z{t),  teS;  let  denote  the  minimum-mean-square  estimate 
of  x^,  given  z(t),  and  22  denote  the  minimum-mean- square  estimate  of  X2,  given  z{t).  Let  x|; 

z'(t),  reS  be  Gaussian  random  variables  with  the  same  covariance  matrices  as  X2  and 
z(t),  tcS.  Then 

ax^'T'bx^  =  ax^”T'bx^  =  E  {ax^  +  bx^|z(T),  reS} 

=  aE{x||z(T),  TfS}  +  bE{x^|z(T)  reS) 

=  axl  +  b$l  =  ax.  +  bx, 
idle 

statement  2;—  Let  x  be  a  random  variable  and  let  F  be  some  class  of  functionals  on  the 
given  data.  Let  2  be  the  functional  of  this  class  such  that 

£  =  E  {(x  —  x)^}  is  a  minimum. 

Let  f  be  an  arbitrary  functional  in  F.  Consider  the  estimate  2  +  ef: 

£‘  +  6f  =  E{(x-2  -  €f)^} 

=  E{(x-2)^} -2eE{f(x-2)}  +  E{f^}  , 

a  f  = -2€  E  {f(x  -  2)}  +  E  {f^l 
The  equation  determining  x  is 

ac  le=0  ° 

To  show  that  x  is  uniquely  determined  by  Eq.  (A-1)  and  truly  yields  a  minimum,  note  that  for 
X  satisfying  (A-1) 

6£=  E  {f^}  ^0  for  f  3fc  0 
Statement  3:—  Define 

Q(t)  =  *(t)  S(t  -  1)  ♦'^(t)  +  W(t  -  1) 

Then  Eqs.  (70)  and  (67)  combine  to  give 

S(t)  =  Q(t)  -Q(t)  H'^(t)[V(t)  +  H(t)  Q(t)  H’^(t)]'^  H(t)  Q(t) 

Using  the  matrix  identity* 

(A  +  BD'S"^)'^  =  A"^  -  A"^B(D  +  b'^A'S)"*  b'^A'^  , 


*  This  identity  can  be  derived  by  using  the  formuios  for  the  inverse  of  a  partitioned  matrix.  See  Frobenius*  re¬ 
lation  in  Ref.  7. 
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we  obtain 


S"^t)  =  Q‘^(t)  +  V'^t)  H(t)  , 

from  which  Eq.  (77)  follows  immediately. 

Similarly,  for  the  predictor  itself,  Eq.  (59)  can  be  written 

£(t)  =  Pj(t)  *{t)  2(t  -  1)  +  C(t)  z(t)  , 


where 

Pl(t)  =  I-Q(t)  H'^(t)  [V(t)  +  H(t)  Q(t)  H'^(t)l”^  H(t)  =  S(t)  Q’^t) 


Now,  byEq.(70), 

C(t)  H(t)  Q(t)  =  Q(t)- S(t)  , 
or,  after  manipulation, 

S‘^t)C(t)  H(t)  =  S‘^t)-Q(t)‘‘  =H'^(t)V'^t)  H(t) 
C(t)  =  2(t)  H'^(t)  V'^t) 


Thus, 

5(t)  =  2(t)  (Q'^t)  ♦(t)i(t-  1)  +  H'’’(t)  V^t)  z(t)]  . 

which  is  Eq.  (78). 
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